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Abstract. This work reports on a study of the spatially coarse-grained velocity dispersion in cosmological N-body simula- 
tions (OCDM and ACDM models) as a function of time (redshifts z = 0-4) and of the coarsening length (0.6-20 /i^^Mpc). 
The main result is the discovery of a polytropic relationship T\ cx (P'^^ between the velocity-dispersion kinetic energy 
density of the coarsening cells, 1\, and their mass density, q. The exponent r;, dependent on time and coarsening scale, 
is a compact measure of the deviations from the naive virial prediction r^viriai = 0. This relationship supports the "poly- 
tropic assumption" which has been employed in theoretical models for the growth of cosmological structure by gravitational 
instability. 
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1. Introduction 

Some recent works 

[Buc hert & Dominguez 1998) 
pommguez et al. "19991 
iMaartens et al. 1999 
iMorita & Tatekawa 20011 



(jBarbero et al. 1997 , 
IMler & Buchert 1999. 
IBuchert et al. 19991 
[Dom inguez 2000 
jPom mguez 2002, 
ITatekawa et al. 20021 have explored a formulation of 
hydrodynamic kind for the formation of cosmological 
structures by gravitational instability. It intends to describe 
the dynamical evolution of the few most relevant fields 
(typically the coarse-grained mass density and velocity 
fields) in terms of a set of autonomous equations for those 
fields, much in the same way as the hydrodynamic equations 
for usual fluids. The widely used dust model ( iPeebles 19801 
ISahni & Coles 1995^ belongs to this class, but it has short- 
comings, most noticeably the emergence of singularities, 
beyond which its application is invalid. An interesting result 
of the systematic study of the hydrodynamic formulation 
are the corrections to the dust model. They arise from the 
nonlinear coupling of the evolution to the structure below 
the coarsening length, and turn out to be relevant in the 
nonlinear regime (^ the regime of large density fluctuations 
about homogeneity). The corrected equations are general- 
izations of the adhesion model (.Kofman & Shandarin 19881 
IGurbatov et al. 1989llKo&ian et al. 1992llMelott et al. 19941 
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[Sathyaprakash et al. 1995|f , which is able to reproduce 
successfully the gross features of the structure evolved by 
gravitational instability in different cosmological scenarios. 
The hydrodynamic formulation requires the corrections 
to dust to be expressed as functions of the coarse-grained 
mass density and velocity fields, so as to get a closed set 
of equations for these fields. This is, however, a major 
theoretical problem, since one cannot follow the standard 
procedure to close the hydrodynamic hierarchy by invoking 
"local thermal equilibrium" (Huang 1987, Balescu 1991D : 
the ideas and the formalism of thermodynamics cannot be 
applied straightforwardly to a system dominated by its own 
gravity, this being in fact still an open question (see, e.g., the 
concise review by Hut JHut 1997t ). 

In this work I report on the empirical search for closure 
relationships with the help of N-body simulations of large- 
scale structure formation. 1 consider in particular the velocity 
dispersion of the particles contained in any coarsening cell, 
which shows up in the equation for momentum conservation 
(in the hydrodynamic parlance, the kinetic contribution to the 
pressure and to the viscous stresses). The simulation results 
for the trace of the velocity dispersion (the internal kinetic 
energy density of the coarsening cells) can be fitted quite 
well by a "polytropic relationship", borrowing the terminol- 
ogy from thermodynamics: X\ cx p'^^''. This result supports 
the polytropic approximation used in the theoretical deriva- 
tion of adhesion-like models cBuchert et al. 1999j . The poly- 
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Fig. 1. Typical log-log scatter plot of the trace of the velocity 
dispersion tensor (in units of km^ s^^ (Mpc//i)^'^, setting 
the particle mass m = 1) vs the number of particles. The 
data correspond to the ACDM model at redshift z = with a 
smoothing length L = 4.55 /i^^Mpc. 

tropic relationship improves with time, and so I conclude that 
it must be a consequence of the evolution by gravitational in- 
stability. For small coarsening lengths or large times, the ex- 
ponent Tj is close to, but significantly different from the virial 
prediction following from the simple model of isolated, re- 
laxed, structureless halos. 

The degree of anisotropy of the velocity dispersion (the 
departure from a spherically symmetric distribution of the 
dispersion) has also been studied, but in this case the scatter 
of the data is large and the results do not provide a significant 
conclusion. The anisotropy tends to decrease with increasing 
mass density (till a few percents at the largest probed densi- 
ties), being however always larger than that associated to a 
Maxwellian distribution with the same density. This decrease 
at least speaks in favor of the approximation of isotropic 
velocity dispersion also employed in deriving adhesion-like 
models. 

This work is organized as follows: in Sec.|2]l detail the 
method employed in the analysis of the simulations. Sec. |3l 
presents the results, which are discussed with the help of the- 
oretical arguments in Sec.|3 Finally, a brief mathematical dis- 
cussion of some topics required in the main text are collected 
in the Appendix. 

2. Method 

The analyzed CDM simulations were performed by The Hy- 
dra Consortium JCouchman et al . 1995 1 using the AP^M al- 
gorithm. They consist of a cubic box (periodic boundary con- 
ditions) of side length 100 /i^^Mpc at the present epoch, con- 



taining N = 86"^ particles. Two different cosmological mod- 
els have been considered, OCDM (fim — 0.3, Q.a — 0.0, 
h = 0.81, CTg = 1.06, r = 0.25, m = 1.63 x 10" Mq) 
and ACDM (f7„ = 0.3, VL^ = 0.7, h ^ 0.96, erg = 1-22, 
r = 0.25, m = 1.37 x 10^^ Mq), at three times, correspond- 
ing to redshifts z — 0, 1.4, 3.6. 

The hydrodynamic formulation can be derived by means 
of a coarsening procedure (Dommguez 2000, 2002). Given a 
comoving smoothing scale L, the coarse-grained mass den- 
sity, velocity and velocity dispersion fields are defined re- 
spectively as follows: 

N . . 

N 

n(x, t) = ~ "(^' ^)] ® 

(n is a second-rank tensor, (g) denoting a dyadic product). 
Here, a{t) denotes the expansion factor, m is the mass of a 
particle, and represent the comoving position and pe- 
culiar velocity, respectively, of the a-th particle, and VF(-) is a 
(normalized) smoothing window. Exact dynamical equations 
can be derived for g and u, expressing mass and momentum 
conservation; the latter equation features the velocity disper- 
sion in the form of a term V • H. The purpose of the present 
study is to check if there exists any approximate relationship 
between 11 and the field g, so that an autonomous set of equa- 
tions can be written for the fields g and u. 

Starting from the coordinates {xq, Ua}a=i...N provided 
by the simulation, the definitions Q were implemented with 
a cubic top-hat window, 

Wiz) = 9{1 - 2|zi|) 6(1 ~ 2|z2|) 0(1 - 2|z3|), (2) 

where 9{-) is the step function. In total 13 different values of 
L were explored, spanning the range from 0.6 h~^Mpc up 
to 20 /i^^Mpc and equally separated in a logarithmic scale. 
For each value of L, at least 2 x 10'' randomly centered, non 
empty coarsening cells were probed. 

The velocity dispersion tensor Ilij{'x) was studied as a 
function of the density g{x.), or equivalently, of the num- 
ber of particles contained in the cell at x, namely n(x) = 
[{aL)^ /m]g{x.). For the purposes of this work, it sufficed to 
consider the three eigenvalues Xi of the tensor 11, or better, its 
three principal scalar invariants: 

Xi = trn = Ai + A2 + A3, 

X2 - ^ [(tr n)2 - tr (n : H)] = A1A2 + A2A3 + A3A1, 

X3 = detn = AiA2A3. (3) 

The quantity (l/2)Xi is the peculiar kinetic energy per unit 
volume due to the motion of the particles relative to the cell 
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3 RESULTS 




Fig. 2. Binned log-log plot of the trace of the velocity dispersion tensor (same units as Fig.0 vs the number of particles 
for the ACDM model. The error bars correspond to the estimated l-cr variance. On the left, for three redshifts (from top to 
bottom, z — 3.6, 1.4, 0) at a fixed smoothing length L = 4.55 /i^^Mpc. On the right, for three smoothing scales (from top to 
bottom, L — 1.41, 4.55, 8.33 /i^^Mpc) at a fixed redshift z = 1.4. The solid line corresponds to the fit (|6j, and the dashed 
line marks Uc, Eq. (|5j- 



center of mass. The other two invariants can be related to 
the degree of anisotropy of the tensor H. More precisely, the 
dimensionless coefficients 

a = l-^, = (4) 

quantify the departure from the isotropic case, Ai = A2 = 
A3 = p. I show in the Appendix that < a, /? < 1, and they 
vanish if and only if 11 is isotropic. Moreover, it follows from 
its definitions that /3 « 3a for small anisotropy. 

To check that this smoothing algorithm was right, it was 
applied to an ideal gas simulation. The results for the depen- 
dence of the kinetic pressure p = (l/3)Ii and the anisotropy 
parameters a, /3 on n agree with the predicted relationships 
(see the Appendix). 

In the next Section I also discuss how the results depend 
on the smoothing details (e.g., on the choice of the window 
(|2jl). I anticipate that the conclusions to be extracted are ro- 
bust. 

3. Results 

The results concerning the three invariants as a function of the 
particle number are qualitatively the same in the whole range 
of smoothing scales and times explored, and for the two cos- 
mological models OCDM and ACDM. Fig.^is a typical ex- 
ample of a Ti vs n scatter plot. The data suggest a polytropic 
relationship between the scalar invariants and n, more and 
more acurate for larger n. Hence, it appears that the increas- 
ingly larger scatter for smaller particle numbers is mainly 



due to the discrete character of the variable n. To eliminate 
this noise, the log-log plots of the raw data were binned into 
40 subintervals on the n-axis: this effectively means averag- 
ing the scalar invariants for different fixed values of n and 
provides also an estimate of the variance. This scatter is in 
fact the main source of error in the parameters 77 and k of 
the fit (|6j. Fig. |2] collects a set of representative cases af- 
ter implementation of this procedure, where the advocated 
polytropic dependence is evident. To be sure that this averag- 
ing method does not introduce artificial features, the analysis 
was repeated by varying the number of bins, from 10 up to 
the limiting case in which each possible value of the discrete 
variable n is treated as a bin. The conclusions are robust and 
the results do not depend on the amount of binning. (Changes 
of the number of bins within this range induce variations in 
the best-fit parameters rj and k which are well within the error 
region in Fig.|3l. The choice of 40 bins is a good compromise 
that efficiently discards the noise but still preserves the rele- 
vant features of the Ti-n relationship. 

Another check probed the influence of the choice of 
smoothing window. The differences found between different 
windows could be explained as a consequence of the dis- 
crete nature of n. Fig. |3l shows the typical result after using 
three different smoothing windows: a cubic top-hat, Eq. a 
spherical top-hat, W{z) = 9{1 - (47r/3)^/^|z|), and a Gaus- 
sian, W{z) ~ exp(— 7r|zp). The windows are normalized 
to unity and give the same weight (=1) to the origin: the 
comparison is then straightforward. The cubic and spherical 
top-hat windows yield indistinguishable results, demonstrat- 
ing that window anisotropy is not relevant (this could also be 
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Fig. 3. Log-log plot of the trace of the velocity dispersion 
tensor vs the number of particles, computed with three differ- 
ent smoothing windows. The data correpond to the ACDM 
model, z = Q,L = 3.45 /i-^Mpc. 



expected, see the last paragraph in this Sec). The Gaussian 
smoothing also gives similar results, but now a slight differ- 
ence in the small-rt end is observed, because the discreteness 
of n implies (i) that 11 = if n = 1 for the top-hat windows 
and (ii) that for small n, 11 is dominated by far contributions 
from the Gaussian tails. However, as discussed later, the dis- 
creteness constraint 11 = if n = 1 for a top-hat window 
can be taken into account in a simple manner and indeed the 
poly tropic fit (|6j holds very well even for n ^ 1. Hence, the 
cubic top-hat window was selected, because it is also most 
easily implemented. 

A final check, carried out only for the particular time z = 
0, was to restrict the coarsening procedure to a subvolume of 
the simulation box (1/64 of the total volume), so as to find 
out to what extent the results are exclusive of the simulation 
volume of 100 /i^^Mpc side length: as expected, the quality 
of the fits are slightly worse because of the reduced number 
of particles, but the conclusions remain the same. 

As exemplified by Fig. |2] one finds a well defined rela- 
tionship between Ii and n, in which three different behav- 
iors can be identified: a polytropic dependence Ii oc ■n?~'^' 
for n larger than a certain value ric, a bending upwards for 
1 < n < ric, and finally a bending downwards for n close 
to 1 which is due to the constraint II = if n = L The in- 
termediate behavior is absent at late times/small coarsening 
lengths, when it is masked by the latter constraint. The value 
of ric increases with the redshift and the smoothing length. A 
rough estimate by eye of the function nc(£, z) is given by 

nc^F{z)nL, ni := (86i/100)^ (5) 

ul is the average particle number in a cell of comoving side 
length L (in units of /i^^Mpc), and F{z) is a mild function 



of the redshift, F ^ 1 for z = 3.6 and decreases slowly 
for smaller z {F = 0.25 for z = 0). This fit for nc{L,z) 
is conservative in the sense of sUghtly overestimating Uc as 
L decreases, however it is precise enough: Fig.|2lshows that 
the polytropic fit is in fact still well followed by the data for 
a certain range below the chosen ric- The best-fit parameters 
7] and K in Eq. (|6j are quite insensitive to the precise value 
of Uc- The intermediate range 1 < n < Uc corresponds 
to underdense cells {ric is associated to a density contrast 
6c = F{z) — 1); the deviations from the polytropic fit may 
then be a remnant of the artificial discreteness, most evident 
in the lattice structure of the initial conditions. Nonetheless, 
this failure of the polytropic relationship is likely unimpor- 
tant from a dynamical point of view, being restricted to ever 
more rarefied regions {5c ~ —0.75 at z — 0). 

The binned data in the range n > Uc were fitted to the 
polytropic relationship 



(aL) 
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(6) 



The factor n — 1 enforces the discreteness constraint Ii = 
for n = 1. This improves the fit on the small-7i region in the 
cases that it extends down to n ~ 1 and it is also suggested 
by the simulated ideal gas: the ideal gas pressure-density re- 
lation is still obeyed with the replacement n ^ n — 1 when 
the smoothing length is so small that most of the coarsening 
cells contain just a few particles. The parameters rj and k were 
determined by the least-squares method. The fit was carried 
out only when the fit range in n spanned at least a decade (in 
which case the range in Ti always extends over at least two 
decades); the polytropic fit is still consistent when this condi- 
tion is not met, but the large errors deprive it of significance. 
In the best cases, the fit range in n included three decades. 
The time and smoothing-length dependences of the parame- 
ters rj and k for the ACDM model are plotted in Fig.|3 The 
results for the OCDM model are very close. 

The other two invariants, T2 and X3, are also fitted by a 
polytropic dependence. In fact, it turns out that the anisotropy 
parameters a and (3 computed with the fitting functions are 
extremely small and 71-independent. The reason is that the 
binned (=averaged) velocity dispersion, being only a function 
of the scalar n, must be isotropic. The anisotropy parameters 
must be calculated for the raw data, before binning: Fig. [S] 
shows a typical a vs n scatter plot. The plots of the parame- 
ter (3 look similar (in fact, the relation /3 « 3a is a very good 
approximation akeady for a ^ 0.1: then, one should rather 
study, e.g., j3 — 2>aX.o gain non-redundant information). Apart 
from the large scatter for small particle number due to the 
discreteness of n, there is also an important dispersion even 
for large ?i, because the anisotropy must be determined by 
something else than only a scalar. Moreover, unlike the "ex- 
tensive" scalar invariants (= sums of positive contributions 
from many particles), the anisotropy parameters, being essen- 
tially the difference between eigenvalues, Eq. \A2\ . are more 
sensitive to the discreteness of n. Thus, the data quality just 
allows to confirm that a and (3 tend to decay with increas- 
ing n but to remain somewhat larger than the average ideal 
gas anisotropy, Ofidcai ~ 5/ (3n) (see the Appendix). No sig- 
nificant dependence (if any) with time and smoothing length 
could be detected. 
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Fig. 4. Plot of the best-fit parameters in Eq. ^ vs the smoothing length for the ACDM model at the three times studied; k in 
units of km^ s~^. The scale of nonlinearity is vq, Eq. 0. The two lines in the 77-plot delimit the estimated 1-cr error region at 
time z = (the errors are of the same size at the other times and they have been omitted for clarity). The error bars for k are 
about the same size as the plot symbols. The virial prediction reads rj — and k cx 1/L (solid line in the K-plot). 



4. Discussion and conclusion 

The simulations show a clear evidence for a polytropic rela- 
tionship (|6j between the density and the scalar invariants of 
the velocity dispersion tensor. The fits improve with decreas- 
ing redshift or smoothing scale; hence one is lead to view 
them as a general consequence of the evolution by gravita- 
tional instability. In agreement with this interpretation, this 
relationship occurs for both the ACDM and the OCDM mod- 
els. Thus, it seems that the existence of the polytropic depen- 
dence itself is independent of the background cosmological 
model, which would perhaps only affect the values of the fit- 
ting parameters slightly. 

The theoretical explanation of the precise values of the 
fitting parameters is not evident. Fig. |3 shows that the func- 
tion 7/(L, z) can be approximately written in fact as a func- 
tion of the single variable L/ro(z), where ro(z) is the scale 
of nonlinearity, defined by the condition that the variance of 
the density contrast is unity: 

^^"-^^'^ ^1, forL = ro. (7) 

The ACDM simulation provides the values tq = 
3.2, 7.7, 16.7 h~^Mpc, respectively for the three considered 
redshifts. This property agrees nicely with the evolution by 
gravitational instability in a hierarchical scenario. The points 
for k{L, z) can be made to collapse on a single curve too if 
K is also suitably rescaled by a z-dependent factor; but this 
must be viewed as pure phenomenology, since no theoretical 
explanation for the values of this rescaling factor could be 
given. 



A theoretical argumentation by Buchert and Dommguez 
( |Buchert & Dominguez 1998 1 provides a polytropic relation- 
ship with a fixed 2 — rj = 5/3: this is equivalent to the 
adiabatic evolution of an ideal gas and was justified for 
early times and under restrictive initial conditions. However, 
Fig- S shows that, precisely in the opposite limit of nonlin- 
ear scales/large times (L < tq), rj is close to this "adiabatic 
value", which even seems to be an asymptote. But the errors 
are too large and the probed range of nonlinear lengths too 
narrow to draw a firm conclusion. 

This flattening of ri{L) suggests another possible expla- 
nation of the polytropic relationship. It seems sensible to hy- 
pothesise that, for smoothing lengths L well in the nonlin- 
ear regime, the kinetic energy should be fixed by the virial 
theorem: if the coarsening cells can be idealized as viri- 
alized, structureless halos, then (aL)^Xi should be propor- 
tional to the potential gravitational energy of the coarsening 
cell, which could in turn be estimated as ~ G{ranY / aL. 
This implies immediately a polytropic dependence (|6j with 
?7viriai = and Kvirial « (flL)^^. This idealized model was 
employed to simulate scatter plots as Fig. ^ assume exact 
isotropy (a = /3 = 0) and the same distribution of particle 
numbers n as obtained in the simulations, and compute Ii 
from the condition that the differences — u in the defini- 
tion Q of n are independent, Gaussian distributed random 
variables with zero mean and a variance given by the virial 
theorem. In this way, an estimated cr < 2 ■ 10^^ for r^viiiai 
was gained. Fig. |3 shows that the best-fit values of rj and k 
are close to the virial predictions for small L; but the devia- 
tions of 77 are well above the estimated fluctuations in ryviriai- 



APPENDIX A: THE ANISOTROPY PARAMETERS a AND /3 
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Fig. 5. Typical log-log scatter plot of the anisotropy param- 
eter a vs the particle number. The data correspond to the 
ACDM model, z = 0, L = 4.55 h^^Mpc. The straight Hne 
represents the average anisotropy for a Maxwellian distribu- 
tion. 

The reason for this discrepancy must lie on the assump- 
tions involved in the above reasoning: (i) the potential en- 
ergy should be dominated by the contribution from the struc- 
ture on the scale L and (ii) each coarsening cell should be 
approximately isolated and in a relaxed, stationary state, so 
that the virial theorem holds. The long-range nature of grav- 
ity could perhaps justify assumption (i) even in a hierarchical 
scenario, when the matter distribution is far from smooth in- 
side the coarsening cells, provided it is not too diluted on the 
scale L either But it cannot be excluded that the corrections 
due to the substructure contribute significantly to the discrep- 
ancy. Assumption (ii) can be easily violated: the coarsen- 
ing cells may have a significant interaction with neighboring 
cells, or be part of a larger virialized halo or contain a bunch 
of streaming particles far from any (quasi-)stationary state. 
As evidence, consider Fig. 3 in dKnebe & Miille r 1999> : it 
is equivalent to my Ti-n plots but the points correspond to 
groups determined by a friends-of-friends algorithm, rather 
than by coarsening boxes of fixed side length. What Knebe 
and Miiller identify as unvirialized groups clearly tend to 
yield a positive ry, in agreement with the trend observed in 
Fig. 0] Therefore, a very interesting result is that the viola- 
tions to the "virialized halo" conditions (i-ii) do not destroy 
the polytropic relationship itself predicted by the virial theo- 
rem, but only change the values of the parameters. 

In conclusion, N-body simulations have provided ev- 
idence for a polytropic relationship between the coarse- 
grained mass density and the peculiar kinetic energy. This 
relation can be characterized by an exponent 77 which mea- 
sures the (scale and time dependent) departures from the 
"virialized halo" prediction. It was found that, for smooth- 



ing lengths well in the nonlinear regime, rj is significantly 
larger than zero. There still remains the task of theoretically 
explaining this polytropic dependence. Future work in this 
direction will address scale-invariant cosmological models: 
they provide testbed cases which can be easily controlled in 
simulations and easily analyzed theoretically. 

Appendix A: The anisotropy parameters a and 

In this Appendix I derive some properties of the anisotropy 
parameters defined in Eqs. ©. Since the tensor II is positive- 
definite, then Xi > 0. This implies the bounds < a, (3 < 1. 
In fact, Ii > yields immediately that a, /3 < 1. On the other 
hand, for any triplet of non-negative numbers, the following 
inequalities hold ( |Abramowitz & Stegun 1965) : 

3 

A,^ > A1A2 + A2A3 + A3A1, 

i=l 
1 ^ 

:r5]A, >(AiA2A3)'/', (Al) 

and the equality is satisfied if and only if all three numbers 
are equal. Combining these inequalities with the definitions 
( I3I4> . it is found that a and [3 must be non-negative, and they 
are zero if and only if the tensor II is isotropic. 

Let us write A.^ = ^li + SXi, so that the coefficients SXi 
represent the deviation from isotropy and satisfy SXi — 0. 
Then: 

a = -^Y.(^X,){SX,), (A2) 

1 i<j 

27 

P^3a-^{SXi)iSX2){SX3), 

and in the limit of small anisotropy, \SXi\ <^ Ti, the equality 
p « 3a holds. 

The simulations of the ideal gas provide (aidcaOn ~ 
5/(3n) with a large scatter, for the reason discussed towards 
the end of Sec.|3l((- • •)« refers to the binning procedure de- 
tailed in that Sec). Nevertheless, this result is reliable be- 
cause the 1 /n dependence can be easily explained: II for a 
coarsening cell containing n particles is the sum of n inde- 
pendent random variables. Hence, for not too small n, the 
average (n)„ will be isotropic and extensive, with the scal- 
ing Ti ^ n, while fluctuations around isotropy will scale like 
\dX\ ^ ^fn. Expression (IA2> then yields that a ~ 
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